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Abstract 

SU(A''c) Yang-Mills theory is investigated at finite densities of Nf heavy 
quark flavors. The calculation of the (continuum) quark determinant in the 
large-mass limit is performed by analytic methods and results in an effective 
gluonic action. This action is then subject to a lattice representation of the 
gluon fields and computer simulations. The approach maintains the same 
number of quark degrees of freedom as in the continuum formulation and a 
physical heavy quark limit (to be contrasted with the quenched approximation 
Nf — > 0). The proper scaling towards the continuum limit is manifest. We 
study the partition function for given values of the chemical potential as well 
as the partition function which is projected onto a definite baryon number. 
First numerical results for an SU(2) gauge theory are presented. We briefly 
discuss the breaking of the color-electric string at finite densities and shed light 
onto the origin of the overlap problem inherent in the Glasgow approach. 



PACS: 11.15.Ha, 12.38. Gc 
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1 Introduction 



The next generation of particle accelerators (RHIC, LHC), which will start operating 
at the beginning of the next millennium, will probe the deconfined regime of QCD, 
the theory of strong interactions, and might reveal exotic states of hadronic matter 
which appear under extreme conditions, i.e., temperature and density. Due to a sig- 
nificant increase in computational power in the recent past, numerical simulations 
of lattice QCD have provided insights into the high temperature and zero density 
phase and have predicted a series of interesting phenomena [ij , such as deconfine- 
ment and restoration of chiral symmetry. Unfortunately, an adequate description of 
finite density hadron matter is still lacking due to conceptual problems in setting up 
an appropriate " statistical" measure which can be handled in computer simulations. 

The generic approach to Yang-Mills thermodynamics at finite densities is based on 
the introduction of a non-zero chemical potential. In the case of a SU(2) gauge 
group, the fermion determinant is real and can be included in the probabilistic 
measure. Numerical simulations can be performed by using standard algorithms p[, 
although this numerical approach consumes a lot of computer time due to the non- 
local nature of the action. Recent progress for the case of a SU(2) gauge group can 
be found in 0. In the case of a SU(3) gauge group, the fermion determinant acquires 
imaginary parts for a non- vanishing chemical potential and cannot be considered to 
be part of the probabilistic measure. The most prominent example to circumvent 
this conceptual difficulty is the so-called Glasgow algorithm There, the fermion 
determinant is considered to be part of the correlation function to be calculated, 
and the probabilistic measure of zero-density Yang-Mills theory is used to generate 
the gauge field configurations. However, it turns out that this approach suffers from 
the so-called "overlap" problem implying that for realistic lattice sizes an unrealistic 
number of Monte-Carlo steps is necessary to achieve reliable results 0. 

In order to alleviate this problem, the so-called quenched approximation, i.e., the 
limit Nf — > 0, where Nf is the number of quark fiavors, greatly reduces the numerical 
task for practical calculations. While in the case of real QCD one expects a drastic 
change in the hadron density for a chemical potential /i = fionset ~ rns/S {ttib is the 
baryon mass), the onset value fionset which is extracted from quenched lattice QCD 
seems to be unnaturally small Subsequently, it turned out that the quenched 
approximation, i.e., the limit Nf 0, of lattice QCD does not meet with the naive 
expectation that this limit coincides with the heavy quark limit of real QCD 0. 

The confining properties of Yang-Mills theory in the desired limit where the quark 
mass and the chemical potential are simultaneously made large (heavy quark limit) 
while the density of quarks is kept finite and non-zero were first investigated in ^. 
Using Kogut-Susskind quarks, this limit simplifies the fermion determinant and 
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allows a significant improvement of the statistics [||. A recent breakthrough [|1^ was 
achieved by resorting to the heavy quark limit in addition to the use of the canonical 
ensemble, i.e., the hadron system of fixed density (by contrast to the grand-canonical 
ensemble of fixed chemical potential). The canonical approach involves the grand- 
canonical partition function with imaginary values of the chemical potential as first 



pointed out in |TT[. A great simplification arises in this case from the fact that 



the fermion determinant is real. The numerical analysis of |jTO|] reveals that the 



deconfinement phase transition becomes a cross-over at finite density (see also p). 
The string between static quark breaks yielding a constant heavy quark potential 
at large distances. 

In this paper, we present a new approach to Yang-Mills theory at a finite density of 
heavy quarks. We shall calculate the continuum fermion determinant for arbitrary 
entries for the gluon field in the large mass limit (rather than in the quenched 
approximation Nf 0) using the Schwinger proper-time regularization. The result 
is a gauge invariant action of the gluon fields, depending on an UV-regulator A, 
and is added to the standard action of Yang-Mills theory. The total result can be 
discretized on a lattice with lattice spacing a by standard methods and can be used 
as input for computer simulations. In the critical limit A —>■ oo, a —>■ 0, physical 
quantities become independent of the regularization scheme, and the results for 
observables are independent of the choice of the technique, i.e., lattice regularization 
or Schwinger proper-time regularization of the fermion determinant. The advantages 
of our approach are as follows: firstly, the starting point of the calculation is the 
continuum quark determinant with the correct number of degrees of freedom. The 
approach is not plagued by spurious states, and its heavy mass limit is manifestly 
the correct QCD limit. Secondly, despite the fact that the total gluonic action is still 
non-local, it is simple enough to allow for fast computer simulations. Thirdly, the 
correct scaling of physical quantities towards the continuum limit is manifestly the 
same as the one proposed by continuum Yang-Mills theory with quarks included. 

The paper is organized as follows: in the next section, we contrast the heavy quark 
limit with the quenched approximation, and discuss the difference in the renor- 
malization group scaling in either case. The calculation of the (continuum) quark 
determinant in the large-mass limit is presented in section |^. We address the renor- 
malization of the coupled gluon quark system and discuss the approach of the con- 
tinuum limit in lattice simulations of the joint system. At the end of section 3, we 
calculate the canonical partition function describing a system with definite baryon 
number. First numerical results for the case of an SU(2) gauge group are shown 
in section 4. The dependence of the quark density on the chemical potential for 
several values of the temperature below and above the deconfinement temperature 
is discussed. The conclusions are left to the final section. 
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2 Heavy Fermions on the lattice 



2.1 Setup 

Our aim is to address the Yang-Mills theory at finite densities of baryons. For this 
purpose, the grand-canonical partition function in Euclidean space time, i.e., 

= J VqVqVA^ exp |- J d^x {Lq + Lym)^ , (1) 

Lym = -^F^M]F;AA] (2) 
Lq = q{x) ip - m - Y\ - ^/^t" , (3) 

serves as a convenient starting point. Thereby, is the usual field strength 

tensor of the gluon field A^, (7 is the Yang-Mills gauge coupling strength, is the 
chemical potential and m is the quark mass. We will assume iVj quark fiavors 
which are degenerate in mass. The convention for the 7-matrices can be found 
in appendix A. A complete gauge fixing is understood in (P with the gauge fixing 
terms included in the measure 'DA^J_. Below, we will employ the lattice version of the 
formulation ([^-0) implying that we need not address the details of the gauge fixing. 
An alternative description of the grand-canonical partition function is obtained by 
integrating out the quark fields i.e., 

= y VA^ Det - m - 7^^^ - i/i7°] exp ^ d^x LyA/| , (4) 

= 1 VA^, exp|-iVySe//[A](A,m,/i) - ^ rf^x LyA^j , (5) 

where A is an ultra-violet regulator. Our goal will be to calculate S'e//[A](A, m, //) 
for large values of the quark mass m. The result will be a gauge invariant functional 
of the gluon fields A^. The joint action, S^ff + / d^x Lym-, will then be discretized 
on a lattice of spacing a and will be subject of computer simulations. In the critical 
limit, a — >^ 0, A — > 00, physical observables will be independent of the type of 
regularization and will approach the continuum result. 

2.2 Heavy fermion versus quenched limit 

One can think of two limits for specifying the heavy quark approximation, i.e., 

^1/2 2^ ^ A ^ m (quenched limit) (6) 

^1/2^ 2^ <^ m ^ A (heavy quark limit) , (7) 

where T is temperature. The string tension a serves in this case as the typical 
energy scale of the pure Yang-Mills system. In the case of the quenched limit (Bf), 
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the heavy quarks are decoupled from the Yang-Mills system. One-loop perturbation 
theory (see e.g. |T2|) then tells us that the continuum limit is approached via the 
scaling 

a\P) = const, exp | - ^ /? | , (8) 

where N^. is the number of colors, (3 = 2Nc/g'^, and M is an arbitrary physical 
quantity of energy dimension one. 

By contrast, if we would like to associate the impact of charm, bottom and top quark 
on the gluonic sector with heavy quark physics, equation must be considered 
as the physical heavy quark limit. In this case, these quark degrees of freedom 
contribute to the critical behavior, and one finds 

aHl3) = con=t. exp | - _ ] - (9) 

A lattice Monte-Carlo simulation of SU(Ai'c) gauge theory with Nf quark flavors 
must recover the scaling towards the continuum limit. In particular, this scaling 
must be obeyed in the physical heavy quark limit (0). 

Note that the quenched limit (P) is formally recovered from (|^) by taking the limit 
Nf 0. The dependence of the bare coupling constant g on the ultra-violet regu- 
lator a provided by must not be confused with the behavior of the renormalized 
coupling (7_r(s) on the renormalization point s, which for instance enters a renor- 
malization flow analysis first proposed by Wilson ||13[. In the latter case and for an 
energy scale s < rric, where rric is the charm quark mass, the "running" of gR{s) 



with s is dictated by the three light, so-called active, quark flavors 



3 Heavy fermions' action 

3.1 The fermion determinant 

The goal of this section is to calculate the fermion determinant 

Dp := Det^^ (z^ - «/i7° - ^ ) , D^. ■= - A^ix) , (10) 

resorting to a 1/m expansion where m is the fermion mass. We will here only study 
the case where the masses of the quark flavors are equal. Nf is the number of 
quark flavors and an UV-regularization is understood in ([l0|). A^{x) is the SU (iVc) 
gauge field and /i the chemical potential. We will use anti-hermitian 7-matrices 
throughout this paper (see appendix A). 
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The determinant D (|ToD is a Lorentz scalar in four dimensions and therefore invariant 
under a reflection of all its vector entries, i.e., id^ — >• — i(9^, — >• —A^, fj, — >• — /z. 
Exploiting the anti-hermitian property of the 7-matrices, we therefore obtain 



Det^//2 
Det^//2 



n - 







°) (nt + 



n 



m . 



This equation shows the familiar result that the quark determinant is real for van- 
ishing chemical potential. Eq. ( pTD also tells us that one generically expects an 
imaginary part of D for real values of the chemical potential while the determinant 
D is again real for purely imaginary entries of /i (note 7 



ot 
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The functional determinant Dp can be represented as a product of eigenvalues, i.e.. 



«/i7 



nt + 



«/i7 



{reg} 

n A. 



Nf/2 



V^„(x° + 1/T,f) 



(12) 

(13) 
(14) 



It is convenient for technical reasons (see also [|T^]) to remove the chemical potential 
H from the operator by a scale transformation of the spinor i.e.. 



nnt ^'^{x) 
<(x° + i/r,f) 



A„^^(x) , 

[-exp(-/i/T)] <(x°,f)) . 



(15) 
(16) 



Introducing the gauge covariant differential Z)^ '■= + lA^ and the field strength 



P^vi i.G., 



we then find 



(17) 



nnt 



flU • 



Schwinger's proper-time method provides a gauge invariant regularization of func- 
tional determinants. In particular, the contribution of the fermion determinant D 
to the gluonic action becomes 



- InD 



N 



f 



d X lim 



°° dr 



tr K{r]x,y) . 



The so-called heat kernel K{t; x, y) satisfies the equation 

^K{T;x,y) + UU^K{T;x,y) = 



(19) 



(20) 



and the boundary conditions (t = x° 



K{T = 0;x,y) = 5{x - y) , 



Tl Tfi 

K[t; it^ + -,x),{ty + j;,y)] 
lim K{T;x,y) 



-fi/T\n—m 



for 

K{T;x,y) 







x,yeV, (21) 
(22) 
(23) 



where V is the space time volume of the (lattice) universe, T is the temperature and 
n,m G Ij. At the present stage, V is considered to be a cylinder with a periodicity 
of 1/T in time direction and infinite extension in spatial directions. Although the 
source 11 11^ which enters the heat equation (^) is hermitian, the desired imaginary 
parts of ([l9| ) will originate (see below) from the non-trivial topology of the non- 
simp ly connected space-time manifold ||15||. 



In order to derive the systematic expansion in powers of the inverse fermion mass, 
i.e., 1/m, we shall adopt the heat-kernel expansion, and resort to the techniques 
reported in the important paper by Ebert and Reinhardt |T^. This expansion is 
generated by the ansatz 



K{T;x,y) = Ko{T;x,y) H{T;x,y) , 



H{T]x,y) = ^hj{x,y)T^ 

j=0 



By definition KQ{T;x,y) satisfies the equation 



d 

— Ko{T;x,y) + 



m 



KoiT;x,y) = . 



A particular solution to this equation is given by ||T6| 

Ko{T;x,y) 



{ 9 it^-ty + af + {x-y)^' 
exp < —m T — 



(24) 



(25) 



(26) 



4r J ' 

where a is an arbitrary constant. The interacting part, H{t; x, y), of the heat kernel 
satisfies the equation 



' d 1 ^ 



-a F 

2 1^'^ i^'^ 



H{T;x,y) = 



(27) 



where = x^ — y^ + oid^Q. Suppose that H{t] x, y) provides a solution to equation 
(P^). Since the gluon field satisfies periodic boundary conditions, i.e.. 



A^{t + n/T,x) = A^{t,x) 



(28) 
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this solution possesses a discrete translation invariance, i.e. 

H{T;{t^ + n/T,x),{ty + n/T,y)^ = H{T;{t^,x),{ty,y) 
A close inspection of (^) then shows that 

h(t; {tr, + n/T,x), {ty,y 



n integer . (29) 



is also a solution of (pTl) if one chooses a = n/T. Equipped with these prerequisites 
we arrive at the central result of this section: the heat kernel which satisfies the 
boundary conditions (|2^) is given by 

1 



K{T;x,y) 



E( 



exp 



-m T — 



{t^-ty + n/Tf + {x-y)^ 



(30) 



At 



H[t; {t^ + n/T,x),y) 

The so-called diagonal parts of the heat-coefficients, i.e., limx^yhj{x,y), can be 
found in |T6[. In order for ensuring proper boundary conditions imposed by finite 
temperatures, the full functional dependence of hj{x,y) on x,y is needed. The 
calculation of the complete heat-coefficients for j = 0, 1, 2 is one of our goals in the 
present paper. The explicit calculation is shown in appendix B. Below, we will only 
make use of ho{x, y) which is given by 



ho{x,y) 



V exp \ —i ^fii^') dx'^ 

I Cxy 



(31) 



where C^y is a straight line connecting the points x and y. The result for hi{x,y) 
and h2{x,y) can be found in appendix B. Inserting (^) into (|19D while employing 
the expansion ( P^ generates the 1/m expansion of the fermionic contribution to the 
gluonic action, i.e.. 



n j=0 

1 r°° dT f o 



(32) 



n 



exp < —m T — 



ATH 



d'^x tr hj {(tx + n/T, x) ,x) . 



For the later discussion, it is convenient to single out the temperature independent 
part n = 0. Performing the r-integration we finally obtain 



^/ = ^ E {m'-'^ r - 2, ^) / d^x trh.ix, 



XI 



m 



4-2j T. 
J- 1 



n 



1 2 ' 

2 m m 



4T2 ' A2 



— d xtihj {{tx + n/T, 



XI ,x) 



(33) 
(34) 
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where the trace extends over color only, and where r(j, x) is the incomplete gamma 
function, 

roc 

T{j,x) = / dss^-^e-' , (35) 

J X 

and where 

hix^y) = ^^-p{-^-f} • (36) 

Since the limit lim^^^o ^ji^, u) does exist, only the term n = ( pBf ) picks up a diver- 
gence. This observation reflects the familiar fact that only temperature independent 
terms are affected by ultra-violet divergences. 



3.2 Renormalization 

In this subsection, we will study the UV-divergence which emerge in (|33D . Since in 
this equation only the diagonal part of the heat coefficients are involved, one might 
resort to the derivation of the diagonal parts presented by Ebert and Reinhardt 
in 



T6|. One finds (see also appendix B) 



tr/io = tr/ii = 0, tr/i2(x,a;) = ^trF^^^F^^. (37) 



The only non-trivial gluon field dependence is therefore induced by the hj, J > 2 
terms. Since the limit Ymvx^Q T{j — 2, x) exists for j > 3, only the term of (^) with 
j = 2 develops a singularity. In fact, one obtains for a large cutoff A 

The divergent part of the fermionic action can therefore be calculated analytically, 

Renormalization is accomplished by absorbing the divergences by an appropriate 
choice of the bare coupling strength g, i.e.. 



lim 

A— >oo 



^ f ji„ rpa rpa _i_ ndiv _i_ ndiv 
Ag"^ J ^ f^u-^ fJ-iy '-'gluon "T ^/ 



finite . (40) 



One loop perturbation theory of Yang-Mills theory without quarks yields 

1 iV 11 A^ 

4? " ^T^V^^""'' ^''^ 
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where A^^^ is the number of colors and /i is an arbitrary renormahzation scale. In- 
cluding Nf quark flavors, we make use of ( PD| ) and demand 



1 11 Nf 1 , A 



In— H i-- In— finite . (42) 

The renormahzation group /3-function, 

can be readily calculated from (^). We recover the well known result 

. _ #1 f _ 0(i) . (44) 



87r2 V 6 3 

The only role of the heavy fermion determinant at this level of the large m-expansion 
is to generalize the renormahzation group /^-function of the pure Yang-Mills theory 
to the situation with Nf quark flavors included. 

We are now in the position to answer the question which raises from subsection 2.2: 
how can the correct scaling towards the continuum limit be observed in lattice 
Yang-Mills theory with Nf quarks included? If we neglect for the moment the finite- 
temperature corrections (0), the inclusion of fermions only affects the parameter 
P — 2Nc/ in front of the Wilson action, i.e., 

Nf N 

/5-/?F = /5-^Y ^^^"^^ ' ^^^^ 
where we have used the freedom in choosing the (fermionic) cutoff for defining 

A2 1 

(46) 
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where a is the lattice spacing, a is the string tension of pure, i.e., Nf = 0, Yang- 
Mills theory and serves as reference scale. Let us assume that we calculate some 
physical mass M in units of the lattice spacing as function of the only parameter 
Since the numerical simulation is carried out with the standard Wilson action 
(with a coefficient changed from (3 to Pf), we recover the familiar scaling behavior 
at large values of P which is at one loop level 

M' a\Pp) = K exp j-jY^ \ , (47) 
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where k is a numerical constant which must be "measured" by the numerical simu- 
lation. Inserting f3p (^) in ( ^7D elementary manipulations of this equation finally 
yield 



M^a^ = K { — -\ expl p } . (48) 

The crucial observation is that o? exponentially decreases with f3 while the slope 
precisely meets with the expectations for a SU(Ai'c) gauge theory with Nf quark 
fiavors. The result of this subsection is also of practical importance: performing 
numerical simulations with f3p as coefficient of the Wilson action and observing the 
standard scaling (|^) automatically implies the correct approach of the continuum 
limit with Nf fermions included. 

Let us assume that we have calculated two physical observables Mi and M2 and that 



we have obtained the proper /3-dependence ( ^8]) for either quantity at asymptotic 



values of /?. The ratio of both observables then becomes 

Ml ~ ^ \K2Ml) 

Using e.g. M2 as reference scale one expresses Mi in terms of the "measured" quan- 
tities Ki and K2, i.e., 

Ml = — M2 . 

/«2 

We finally note that temperature dependent terms (0), which we have not consid- 
ered in this subsection, are finite in the limit a — > 0. We therefore do not expect 
that these terms will change the critical behavior for (3 —>■ 00 in (^8|). 



3.3 Lattice YM-theory with finite chemical potential 

In this subsection, we will consider the temperature dependent part of the sum (P^. 
The pre-factors of the heat coefficients hj{{tx + n/t, x),x) are given by the functions 
/j(n^m^/4T^, m^/A^) which stay finite when the UV regulator is removed (A 
00). Using the heavy quark, but non-quenched limit m? / (see discussion in 
subsection 2.2), one finds 

2m^T^ / nm\ 

{2mTy \Y) ' 

where the Kn{x) are the modified Bessel functions of the second kind. The heat 
coefficients hj are functionals of the gluonic fields only and carry energy dimension 
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2j. One therefore expects that their order of magnitude is given by Ay^ where 
Aym is the characteristic energy scale of pure Yang-Mills theory. This implies that 
the sum over j in (|^) generates the desired heavy quark expansion in powers of 
Ay^j/mT. In the following, we will assume that the mass is large enough for a 
truncation of this sum at j = 0. 
Using (^) one might interpret 



/io((tx. + l/T,f)) 



P(x) 



(50) 



as Polyakov line P{x) which starts at the space time point x and which winds around 
the torus in time direction finally ending at the same point x. The operator P{x) can 
be directly translated to the corresponding lattice operator: it is the (path-ordered) 
product of the link variables U^{x) along the = 4 direction. Expressing the field 
strength squared in terms of the plaquette variable P^u{x), i.e., 



— tr P,,Jx] 



(51) 



we finally obtain the lattice action Siatt of SU(2) Yang-Mills theory with Nf heavy 
quark flavors 



{^l>u){x} 

{x} ^ " n=l 



kJ—\ 



\T J 

^rP"(x) + e-"T tr {p\xY 



(52) 
(53) 



where I3p is given by (|45|). Exploiting the heavy-quark limit, i.e., m T, one uses 
the asymptotic expression for Bessel functions. 



Koix] 







1 H 









(54) 



In this limit —Sp in (|53|) becomes 



E 

{x} 



-ir 



3/2 



rp OO 

— V 

2m ^5/2 



gnii^ trP"(x) + e"""^ tr (P^(x)") 



(55) 

For a SU(2) gauge group, this equation is real even at finite values of the chemical 
potential as it is trP. By contrast, one expects imaginary parts for real and 
SU(A'c > 3). Finally, one observes that for purely imaginary entries n = iu, u real, 
equation ( 



is real as expected (see discussion in subsection |3.1|) 
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For fi < m, the sum in (p5D is rapidly converging, while the sum is only asymptotic 
for fi ^ m (for a discussion of this issue in the context of QED see [T^). In 



the later case, resorting to different representations of the Bessel function K2 it is 
possible to perform an analytic continuation from the finite sum for fi < m to large 
values of yU > m |T^. The outcome of this procedure is that by means of analytic 



continuation one can assign a finite and real value to the sum also for /i > m. The 
lack of an imaginary part possesses a physical interpretation: the chemical potential 
/i describes the gain in energy if a particle is added to the system, while an energy 
E is necessary to produce the particle. Neglecting binding and confinement effects 
(for this argument only), the particle must occupy an empty phase space cell and 
carries a momentum larger than the Fermi momentum. The loss in energy due to 
particle production is therefore always larger than the gain in energy. The system 
is stable. 

The production of quarks at /i > m is conflicting with quark confinement, which 
does not tolerate single quark productions. This interplay between these mechanisms 
can be anticipated from (|55|) . One expects a drastic influence of the heavy quark 
corrections to the zero density action if 

exp{^} > 1. (56) 

For small temperatures (and moderate densities), {\P\) is small, and the onset of 
the density effects is postponed to values fi ^ m. On the other hand, at high 
temperatures (T > Aym), {\P\) is of order one by temperatures effects only, and 
a significant rise of density is expected for ^ m. From a physical point of view, 
this rise of density might be interpreted as single quark productions which become 
feasible in the deconfinement region. 

In the following, we will study the case where is only slightly larger than m. This 
procedure will allow us to study baryon matter at moderate densities. We hope that 
in this case the truncation of the sum aX n = 1 already captures the essence of the 
asymptotic series (|55|) . The extension of the considerations to larger values of the 
chemical potential is an interesting task which is left to future investigations. In the 
present paper, we will only perform a consistency check by estimating the n = 2 
term. Within this approximation, we finally obtain 

{x} ^ tJ->'^ 

+ — 1= — ■ — a e T trPlx) + e t trP'lx) 
The partition function which describes the interaction of SU(2) gauge fields with Nj 
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heavy fermions is finally expressed as an integral over the link variables U^{x), i.e., 

Zifi) = I VU^{x) expi-Siatt} . (58) 

Equations ( |57| , |58| ) are the main results of the present paper. Note that the consid- 
erations of heavy fermions induce non-local terms to the effective action of gluons. 
These terms are given by the Polyakov line in (p7[). This non-locality is confined 
to the temporal direction while the action is still local in spatial directions. Note 
further that trP(x) is real for a SU(2) gauge group. In the SU(2) case the action in 
(0) can be therefore easily simulated with a moderate increase in computer time 
compared with the case of pure YM theory by resorting to the standard heat-bath 
algorithm. First results of such a simulation will be presented below. 



Let us finally calculate the baryon density p from the functional integral ([58|). The 
derivative of the partition function with respect to the chemical potential yields 

where B is the number of Baryons which are present in the (lattice) universe. In- 
serting (|58| ) and (0) into (|59D, one finds for the baryon density 

P = B/V = j^^^, [e— (trP) - e-— (trPt)] , (60) 

where V = N^a^ is the volume and where T = l/Nta was used. and Nt denote 
the number of lattice point in spatial and time direction, respectively, (tr P) is the 
expectation value of the Polyakov line. For positive values of the chemical potential 
and fi ^ m, thermal excitations of anti-quarks can be neglected. In this case, the 
second term on the right hand side of (|60|) can be neglected, and the density is 
triggered by the expectation value of the Polyakov loop. 



3.4 Lattice YM-theory at fixed baryon density 

Endowed with the results of the previous subsection, it is an easy task to apply 
the approach of Miller and Redlich |TT[] for describing the Yang-Mills theory with 
a definite value of the baryon density. Their pioneering approach is based on the 
observation that the partition function 

Z,{B) = expj-^z.^} Z{zv) (61) 

describes the Yang-Mills theory at a fixed density provided by B baryons in the 
universe. Thereby, Z{iv) is the grand-canonical partition function (P with an imag- 
inary entry as chemical potential. Inserting (|l|) into (|6l|) , one readily verifies that 
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the z/-integration constrains the parameter B to the baryon number, i.e., 

B = ^ Jd'x(^iq{x)j%{x)'j , (62) 



where an average over a time period of length 1/T is understood in (p^). Note that 
T denotes the temperature, and that we use anti-hermitian 7-matrices (see appendix 
A). 



One easily repeats the derivation of the heavy quark lattice action ( p2D for the case 
of an imaginary chemical potential. The calculation of the eigenmodes of the heat 
kernel (subsection is unaffected by the choice of an imaginary chemical potential, 
and one essentially observes a change of the boundary conditions ([T6|) of the quark 
fields. The final result of this consideration is that the heavy quark lattice action 
(0) is given in the case of an imaginary chemical potential by replacing fi iu in 

In accordance with P, one immediately realizes that the heavy quark lattice 
action (|S3D is real for an imaginary chemical potential. In agreement with [^, we 
also observe that the partition function (^) for a given baryon number B is invariant 
under a center transformation of the links which belong to a spatial hypercube at a 
given time slice. In this case, the Polyakov loop acquires a phase exp(i 2tt/Nc) which, 
however, can be absorbed by a redefinition of the integration variable u v—2t{/Nc. 

By contrast to the case of a finite real value of the chemical potential, the sum over n 
in ( |55D is rapidly converging for /i = iz/ as long as m > T. Calculating the fermionic 
contribution exp{— S'^?} to the probabilistic weight from ([55|) , one observes that the 
terms P"(x) with a multiple, i.e., n > 1, winding of the Polyakov loop around the 
torus are suppressed by the factor l/n^^^ in (|55|) . One finds 



e"^^ ~ + C[e*'^''trP(x) + e-'"'^ ii P\x)\] , (63) 

{x} 

Inserting ( |55D into (pT]), the Fourier integration can be explicitly performed. The 
final result is 



(65) 



C^^^ E n trP(x.)l expE(/?F E^.^(^)) 

{x} i=l {x} ^ tJ->y ^ 



where Xi 7^ for z 7^ /c holds. Insertions of P(xa)P^(xb) are suppressed by a factor 
and do not contribute to the leading order result (|65D. From a physical point 
of view, these insertions correspond to particle anti-particle excitations which are 
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suppressed by a factor of roughly exp{— 2m/T} in the probabihstic weight. Note 
that Zp{B) is manifestly center symmetric since only products of Polyakov loops 
which consist of a multiple of Nc factors are present in (|65|) . 

It is instructive to compare the partition function Zp{B) of a system of finite baryon 
density, i.e., B 0, with the zero density result Zp{B = 0). One observes that 
the finite density partition function is suppressed by a factor (^'^"^ relative to the 
zero density case. This illustrates the overlap problem which rises by resorting to 
the Glasgow algorithm Simulating the partition function of a system with B 
baryons employing an algorithm with generates an Monte-Carlo ensemble which is 
based on the zero density partition function requires a huge amount of statistics to 
overcome the entropy factor C^'"^- 



It was argued by Svetitsky ||20[ quite some time ago that the Polyakov line corre- 
sponds to a field configuration with an essential overlap with the wave function of 
a static quark. At that time, this approach paved the way for the interpretation 
of the Polyakov line expectation value as an order parameter for the deconfinement 
phase transition. Our result (|65D nicely confirms Svetitsky's considerations on a 
quantitative level. 

We finally point out that in a recent publication ||2l| a close relation of the phase of 



the (lattice) fermion determinant and the imaginary parts of the Polyakov loop was 
observed. Our results confirm these findings. 



4 Numerical simulations 

While the derivation of the effective gluonic action in the previous subsections holds 
for a SU(A*'c) gauge group, we will below confine ourselves to the case of a SU(2) 
gauge group. The partition function (|58|) supplemented with the action 

{x} ^ I J 

^ = ^^3/2 ^ b " + e- T \ (67) 



can be simulated with the standard heat bath algorithm of Creutz |T9[. The non- 
locality of the action due to the Polyakov loops only yields a modest increase in 
computational time. We adjusted the temperature T of the system by varying the 
number Nt of lattice points in time direction, i.e., T = 1/Nta. In order to remove 
the superficial renormalization point, i.e., j3f, dependence from physical quantities, 
we assume that the one-loop scaling (|47| ) is a good approximation for moderate j3f 
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Figure 1 : The 
number B of Baryons 
which are present in 
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cal potential n in units 
of the heavy quark 
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values, i.e., /?/ G [2.1,2.5]. In order to fix the overall scales, we used 

a{Pf = 2.3) = 0.16 fm , (68) 

which would correspond to a string tension of a = (440 Me V)^ in pure Yang-Mills 
theory (Nj = 0). The density of nuclear matter is roughly given by p ~ 0.15 fm~^. 
This value corresponds on average to 6 x 10~^ Baryons in an elementary cube of 
size and roughly one baryon in a lattice universe of size 12^. 

In practice, it turned out to be convenient to run the simulation for a definite set of 
^-values, e.g. ^ G [0, 1/2], and to calculate the relation between the parameters /x, 
m, T afterwards. The derivation of (|66|) relies on a truncation of the series (^) at 
n = 1. For being consistent the inequality 

e-^^ \,.r ".r (1^1) < 1 • (69) 



25/2 gM/T _^ e^/T 
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must be satisfied. In the present paper, we tacitly assume that the truncation at 
n = 1 is reasonable if ( pD| ) is satisfied. Investigations beyond this approximation are 
left to future work. 

4.1 Quark genesis versus quark confinement 

In subsection |3]^, we have seen that in the deconfined phase ((P) = C(l)) quarks 
can be produced in large numbers if the chemical potential exceeds the value of the 
quark mass. By contrast, in the "confined" phase, i.e., (P) ~ 0, the production 
of single quarks by means of the chemical potential is forbidden and only thermal 
excitations of baryons (i.e., diquarks for the present case of an SU(2) gauge theory) 
are possible. A large number of baryons is expected to occur for n > 2m — b, where 
m is the quark mass and b is the binding energy of the diquark system. In order to 
study this interplay between quark production and confinement we have calculated 
the number B of baryons which are present in the lattice universe as a function of 
the chemical potential for several temperatures. Pure Yang-Mills theory {Nf = 0) 
possesses a second order deconfimement phase transition for T = Tc ^ 300 MeV. 
The result of the simulation using a 12^ x Nt lattice and (3f = 2.3 is shown in 
figure 0. For definiteness, we used m = lOy^. We have checked that the inequality 
(|69|) is satisfied for the parameter ranges producing figure ^ For temperatures 
below Tc, we clearly observe an onset value fionset larger than the quark mass. For 
temperatures larger a significant rise in the baryon number is observed for /i ^ 
m. While in the previous case the net baryon number is produced via baryonic 
excitations, excitations of single quarks contribute to the net baryon number in the 
high temperature phase. 

4.2 String breaking at finite density 

At zero baryon (quark) density, a convenient order parameter O of confinement is 
constructed from the Polyakov line P{x), 




Keeping in mind that the Polyakov line reverses its sign under a center transforma- 
tion, i.e., P{x) — > — P(x), (P(x)) = is mandatory for a realization of center sym- 
metry. In this case, quark confinement is inherent and one observes O oc l/W, 
where V is the space volume. If the temperature exceeds the critical temperature Tc 
(~ 300 MeV for a pure SU(2) gauge theory), a non- vanishing expectation value of 
the Polyakov line, i.e., (P(x)) ^ and therefore O = 0{1), signals the spontaneous 
breakdown of center symmetry and deconfinement. 
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Figure 2: schematic plot: the confinement order parameter (^) of pure SU(2) gauge 
theory in the case of a finite-density system as a function of the baryon density p 
(left panel) and the chemical potential /i (right panel), respectively. 



While at zero density the color electric string between static quark sources might 
extend to arbitrary length, one expects at finite densities a maximum length Zq of 
the color electric string which is controlled by the average distance between the 
background quarks. In this case, string breaking occurs and the heavy quark poten- 
tial saturates at a quark distance Iq. A precise definition of a deconfinement order 
parameter is cumbersome in the finite density case, and has been under debate for 
more than twenty years (for recent progress see p2[). 



Although one does not expect a sharp deconfinement phase transition as a function 
of temperature at finite densities, the gluonic state might undergo drastic changes 
which are triggered by temperature effects and which are signaled by significant 
(but smooth) changes in observables. To these respects, the behavior of O is of 
particular interest also at finite densities. In the present case of an SU(2) gauge 
group (tr P = tr P^) and a system of heavy quarks, the density (|60|) is proportional 
to the expectation value of the Polyakov loop. One therefore finds 

IpI (X KtrP)| < (|trP|) , (71) 

and concludes that O = 0{1) as long as p 7^ 0. Our result nicely confirms the 



observation in []10[ that color electric string breaking occurs as soon as the baryon 



density is non-zero. 

We are led to the following behavior of O as function of the chemical potential and 
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density, respectively: for temperatures significantly below the critical temperature 
and for /i < fionset, the baryon density is practically zero yielding small values of O. If 
the chemical potential /i exceeds fionset, the drastic rise of the density is accompanied 
by a strong increase of O. At large temperatures (T > Tc), O is non-zero due to 
temperatures effects and changes in O due to density effects are moderate in this case 
(see figure Numerous numerical results confirm the qualitative behavior shown 
in figure |^. Quantitative details are not interesting since they strongly depend on 
the actual choice m ^ y/a and on the fine tuning fj, — > m+. 

5 Conclusions 

Our central idea of the present paper is to combine analytic methods for calcu- 
lating the fermion determinant with the lattice description for deriving a valuable 
description of the Yang-Mills system at finite densities. Removing the ultra-violet 
regularization (A — oo in the case of the quark determinant and a — >■ in the case 
of the gluonic functional integral), physical quantities are independent of the type of 
regularization and approach a unique result. We have shown in section 3.2 how the 
proper scaling towards the continuum limit is obtained in the present case of inter- 
est. The advantage of our approach is twofold: firstly, by construction the approach 
is not plagued with spurious quark states (see e.g. 0). Secondly, the physical heavy 
mass limit, i.e., Aym, T <^ m ^ A, is manifest. 

In the case of the grand-canonical partition function, the chemical potential is chosen 
of the order of the quark mass in order to produce significant effects in the quark 
density 0. Our first numerical results for the case of an SU(2) gauge group have 
been presented in section 4. Rising the temperature above the deconfinement critical 
one we observe a decrease of the onset value of the chemical potential at which a 
rapid increase of baryon density is observed. We interpret this result as follows: at 
low temperatures, baryon density is generated by the chemical potential only via 
the production of baryons (i.e., diquarks in the case of an SU(2) gauge group). At 
high temperatures, by contrast, the excitation of single quarks contributing to the 
density becomes feasible due to deconfinement. 

In the case of the canonical partition function Z{B), describing a system with baryon 
number B, we observe that the quark determinant is expressed in terms of products 
of Polyakov loops. In agreement with the findings in |]10|, the determinant is center 
symmetric, and a non-vanishing expectation value of the Polyakov loop at finite 
densities occurs via string breaking ||T0[| . 

Although the calculation of the quark determinant in the present paper is tied to 
the Schwinger proper-time approach which generates the large mass expansion, the 
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basic idea of combining an analytic calculation of the determinant with a subsequent 
lattice representation of the gluon fields is quite universal. Estimating fermion 
determinants by resorting to different types of approximation schemes has a long 
history in the literature. The idea of studying the opposite limit m ^ (chiral 
limit) by applying this new idea seems very appealing to us. 
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A Notation and conventions 

The metric tensor in Minkowski space is 

g^, = diag(l, -1,-1,-1) . (A.l) 
We define Euclidean tensors T(^e) from the tensors in Minkowski space T(m) by 

where r and s are the numbers of zeros within {/xi . . . /iat} and {ui . . . z/„}, respec- 
tively. In particular, we have for the Euclidean time and the Euclidean metric 

^Ie) = ^^Im)^ 9{e) = diag(-l,-l,-l,-l) . (A.3) 

Covariant and contra- variant vectors in Euclidean space differ by an overall sign. For 
a consistent treatment of the symmetries, one is forced to consider the matrices 
as vectors. Therefore, one is naturally led to anti-hermitian Euclidean matrices via 

o, 

7(E) = ^7(M) , 7(E) = 7(M) • (A.4) 

In particular, one finds 

{iQ^ = -l^E). bU^iE)} = ^9iE) = -25... (A.5) 
The so-called Wick rotation is performed by considering the Euclidean tensors ([A.2|) 
as real fields. 

In addition, we define the square of an Euclidean vector field, e.g. V^, by 

:= V^V^ = -V^V . (A.6) 

This implies that V"^ is always a positive quantity (after the wick rotation to Eu- 
clidean space). 
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(a) (b) 

Figure 3: Possibilities (a) and (b) for calculating the holonomy along the triangle 
path shown in (a) (non-Abelian Stokes theorem). 



B The ofF-diagonal heat coefficients 

The aim of this subsection is to calculate the full space time dependence of the heat 
coefficients hk{x,y). We will show that we recover the well known result for the 
diagonal part hk{x,x). For this purpose, our starting point is the recursion relation 

{k + z^D^)hk{x, y) = {D^ + ^-cr^^[D^, Dj\)hk-i{x, y) , (B.l) 

while the equation for hoi^x, y) is given by 

z^Df,ho{x, y) = , := x^ - y^, . (B.2) 

We first show that the gauge covariant connection 



ho{x,y) = U{x,y) = Pexp|-z 
where the path Cxy is a straight line connecting the points x and y, provides a 



A,{x')dx'\ (B.3) 



solution to equation (|B.2|) . For a proof we calculate the holonomy along the triangle 
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path shown in figure ^ By a comparison of the results obtained by using the 
equivalent paths (a) and (b) in figure ^ at the level 0{dx), one finds 

U\x,y)D^U{x,y) = -i dt t T^^{zt + y,y) , (B.4) 

Jo 

T^,{x\y) ■= U\x\y)F^,{x')U{x\y) . (B.5) 

In the Abelian case J^^u{x',y) becomes independent of y and coincides with the 
standard field strength F^y{x'). It is can be easily checked that both sides of ( |B.4D 
transform homogeneously under gauge transformations. Using ( |B.5|) and the anti- 
symmetry of the tensor T^y{x' ,y) under an exchange of the indices and one 
immediately observes that z^D^ho^x.y) = 0. 

In order to calculate the heat coefficients hk{x,y), we decompose 

hkix,y) = U{x,y) gk{x,y) , with 9oix,y) = 1. 
Making extensive use of ( |B.4| ) and the relation 

D^(u{x,y) gkix,y)^ = (^D^U{x,y)^ gkix,y) + U{x,y) (^d^gkix,y)^ , 

one can cast ( p.l|) into a recursion relation for gk{x,y), i.e., 
{k + Zf,df,)gkix,y) 

= ld^-—f z'dx'^Tf,a{x\y)df, 

\ Z J Crr.ll 



Z^ JCx 



2^ JCxy 



^ 2^2 Jc^ 



z'dx'^J^f,a{x',y) / z"dx"f^J^^(3{x" ,y) 
/ z''^dx^&^T^^{x\y) -'^ z'^dx'^T^y{x\y) 

" Cxti Z J Cxy 

z'dx'^ I z"dx'^[J='^a{x',y),J^uf3{x",y)] 

^ Cxy 

/ z''^dx'^{d'^J='^a{x',y) - dlJ^f,a{x',y)) gk-i{x,y) (B.6) 

^ Cxy J 

In particular, the equation for g\ [x, y) becomes 

{1 + z^d^)gi{x,y) = — \f z'dx'^J='f,a{x',y) f z" dx'^J^ ^^{x" , y) 

Z- J Cxy ^ Cxy 

--2 z'^d^'ad'^^t^a{x',y) z'^dx'^Ty,y{x\y) 

Z J Cxv Z J Cxv 



/ z'dx'^ / z" dxl[Tf,a{x' ,y), J^yi3{x" ,y)] 

Cx y Cx y 

f z'^dx'^{d'^r,^{x\y) - dlJ^,Ux',y)) (B.7) 
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Referring to the particular solution of the equation {k + z^d^)A[x, y) = B{x, y) given 

by 

A{x,y) = ^f z'^-^ z'^'B{x,y)dx'i' (B.8) 

we finally obtain 

g,{x,y)^--[ z'-'z',dx'J^ f z"dxlT^{x\y) f z'" dx"^T,p{x'\y) 

Z JCxy \Z JC^ly -^^x'y 

+^2 L z"Hxy;,T,^{x\y) 

z'Ldx"^T,A^\y) (B.9) 
z"dxl J z"'dx"^[T,^{x", y),T.p{x'\ y)] 

2z'^ Jc , y) - K:f^{^\ y))) 

Not showing terms which vanish if the trace over Dirac indices is performed, we find 

h\{x,y) = -V exp I -I f Ai^{x') dx'^\-f z'~h'^ dx\ 

\ [ z" d< J'.aix", y) f z'" dx'H J'.pix'", y) + ... .(B.IO) 

For calculating g2{x,y), we introduce 

x' — zp + y, x" — z't + y — zpt + y, x'" — z's + y — zps + y 
and rewrite gi{x,y) as 

gi{x,y)^- P^dpJ^ tdt sdsZaZi3Ti^aizpt + y,y)Ti^p{zps + y,y) 

tdtZadi^J^i,a{zpt + y,y) 

-'^fif dp tdtTf,„{zpt + y,y) 

[ p^dp [ tdt [ sdszazpi^i^aizpt + y, y),J^^p{zps + y, y)] 
Jo Jo Jo 

dp tdtZa{d^J^„a{zpt + y,y)- d^Tf,a{zpt + y, y)) (B.ll) 



2 
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where d^J^^^{zpt + y,y) = d^J^^^{x,y%=^pt+y Using ( pl6|) , (pISD and (|B.11|) , the 
final result of a lengthy calculation is 

5(2(x, ?/)= ^ /" z'^ dx'^ |- /" p^dp f tdt I sds(2T^y{z'pt + y, y)T^u{z'ps + y) 

+'^z'pd')^{{J^^x{z'pt + y, y),J^^f3{z'ps + y, y)}) 
+ZaZ'^d'^^f,a{z'pt + y, y)J^^p{z'ps + y, y)^ 

+2 (^f dp f tdtT^,{z'pt + y, y) 

+T, I P^dp [ tdt f sdsz'^z'Jj='f,a{z'pt + y, y),Typ{z'ps + y,y)] 
I Jo Jo Jo 

~\fo '^^fo ^'^^^'c.K^-c.iz'pt + y,y)- d,T^^{z'pt + y, y))^ | +(B.12,) 

where we have not shown the traceless terms. In the diagonal limit x ^ y, the 
formulas for gi{x,y) and g2{x,y) greatly simplify and we recover the familiar result 

gi{x,x) = 

= V. (B.13) 

D 
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